% 上图: LE1, LE2, LE3 随 mu 的曲线
% 下图: x 的分岔图    

clc; clear; close all;

mu = 3.45;                            % 固定 mu
k_values = linspace(0.1, 0.46, 1000); % k 扫描范围，采样 1000 个点(可调整)
x0 = 0.1; y0 = 0; z0 = -1;            % 初始条件
nTransient = 1000;                    % 丢弃瞬态步数
nLE = 1000;                           % 用于计算 LE 的步数
nStore = 100;                         % 分岔图中每个 k 存储多少点

% 用于存储 LE
LE1 = zeros(size(k_values));
LE2 = zeros(size(k_values));
LE3 = zeros(size(k_values));

% 用于存储分岔图数据
bif_k = [];
bif_x = [];

% 扫描 k
for i = 1:length(k_values)
    k = k_values(i);

    % (a) 丢弃瞬态
    state = [x0, y0, z0];
    for n = 1:nTransient
        [dx, dy, dz] = mclm(state, mu, k);
        state = [dx, dy, dz];
    end

    % (b) 计算 Lyapunov 指数
    LE_vals = LEs(state, mu, k, nLE);  % 返回 [LE1, LE2, LE3]
    LE1(i) = LE_vals(1);
    LE2(i) = LE_vals(2);
    LE3(i) = LE_vals(3);

    % (c) 分岔图：再做 nStore 次迭代，记录 x
    s = state;
    for mIdx = 1:nStore
        [dx, dy, dz] = mclm(s, mu, k);
        s = [dx, dy, dz];
        bif_k(end+1,1) = k;
        bif_x(end+1,1) = s(1);   % 记录 x
    end
end

% 上方：LE1, LE2, LE3 随 k 的变化
subplot(2,1,1);
plot(k_values, LE1, 'r', 'LineWidth',1.2); hold on;
plot(k_values, LE2, 'm', 'LineWidth',1.2);
plot(k_values, LE3, 'b', 'LineWidth',1.2);
xlabel('k');
ylabel('LE');
legend('LE_1','LE_2','LE_3','Location','best');
title(sprintf('LE Spectra vs k,  \\mu=%.2f, (x_0,y_0,z_0)=(0.1,0,-1)', mu));
grid on;

% 下方：分岔图 (bifurcation diagram)
subplot(2,1,2);
scatter(bif_k, bif_x, 1, 'magenta', 'filled');  % 用绿色点
xlabel('k'); ylabel('x');
title('Bifurcation diagram');
grid on;
